Introduction

Dans le monde de la data science et de l’IA, on parle souvent d’apprentissage automatique, de modèles de boîtes noires ou encore d’interprétabilité. A partir du 20e siècle, les chercheurs ont commencé à élaborer des méthodes permettant de comprendre le comportement et l’interprétation de ces algorithmes fonctionnant de manière opaque. Mais si ça fonctionne, pourquoi l’interprétabilité est-elle si importante ? Nous allons détailler cette notion à travers un exemple.

Supposons qu’une banque utilise un modèle d’apprentissage automatique pour décider d’accorder ou non un prêt à un individu. Ce modèle prend en compte des caractéristiques telles que l’âge, le nombre d’enfants, le statut matrimonial et d’autres variables pour estimer la probabilité qu’un individu rembourse correctement le prêt. Si cette probabilité est de 70 %, comment savoir si cette décision est juste ? Quels facteurs influencent ce pourcentage ? Le statut matrimonial ? Le nombre d’enfants ? Est-ce que certaines caractéristiques sont surévaluées en raison des tendances dans les données d’apprentissage ? Et qu’en est-il des individus dont le cas est atypique ? Peuvent-ils être injustement rejetés en raison de leur singularité ? L’aspect boite noire de l’algorithme soulève une multitude questions.

A travers ce document, nous allons aborder les techniques d’interprétabilité des forêts aléatoires avec les packages randomForestExplainer et pdp.



Les données

Les particules suspendus dans l’air proviennent de diverses origines, aussi bien naturel qu’anthropique. Ces particules ont des compositions chimiques très variable. L’observatoire de la qualité de l’air en Haute-Normandie, Air Normand, mesure la qualité de l’air tout les quart d’heure à travers son réseau d’une douzaine de stations. Ces données sont accesibles via le package VSURF.

Nous allons utiliser les données enregistrés par la station gui.

Description du package VSURF : Le package VSURF est basée sur les forêts aléatoires. + description

suppressMessages(library(VSURF))
# chargement du jeu de données gui
data("gui", package = "VSURF")


Le jeu de données est composé de 18 variables qui sont les suivantes :

  • PM10 : Concentration moyenne journlière de PM10, in μ g/m^3

  • NO, NO2, SO2 : Concentration moyenne journlière de NO, NO2 , SO2, in μ g/m^3

  • T.min, T.max, T.moy : Température minimale, maximale et moyenne journaliers, en degrés Celsius

  • DV.maxvv, DV.dom : Vitesse maximale journaliers et direction du vent dominante, en degrés (pour la direction du vent, 0 degré correspond au Nord)

  • VV.max, VV.moy : Maximum et moyenne journaliers de la vitesse du vent, en m/s

  • PL.som : Pluies journalières en mm

  • HR.min, HR.max, HR.moy : Minimum, maximum et moyenne journaliers de l’humidité relative, en %

  • PA.moy : Pression atmosphérique journalière, en hPa

  • GTrouen, GTlehavre : Gradient de température journalier, en Celsius


Ces données seront utilisés dans le cadre de la prévision de la pollution de l’air par les particules en suspension PM10 en fonction des autres variables météorologiques. La variable cible est PM10.


Aperçu des premières lignes de gui
PM10 NO NO2 SO2 T.min T.max T.moy DV.maxvv DV.dom VV.max
13 42 48 8 -0.6 4.4 1.173913 160 180.0 8
13 22 34 7 -1.3 5.6 2.012500 20 22.5 7
21 44 55 13 -4.2 -0.7 -2.933333 40 180.0 3
DV.dom VV.max VV.moy PL.som HR.min HR.max HR.moy PA.moy GTrouen GTlehavre
180.0 8 5.347826 11 93 97 96.21739 1010.930 0.7 -0.3
22.5 7 4.958333 3 78 96 87.87500 1016.975 -0.5 0.0
180.0 3 2.250000 0 63 91 82.62500 1024.088 2.4 0.1


Traitement des données

summary(gui)
gui <- subset(gui,!is.na(PM10))

Toutes les observations où la variable cible est nulle ont étés supprimés car ces observations n’apporterons aucune informations. Il reste encore 49 observations avec des données manquantes, ce qui représente 4.46% des données d’origine.

Nous allons créer une sous version de gui sans données manquantes que l’on appelle gui_without_na.

gui_without_na <- na.omit(gui)


Analyse rapide des données

Regardons la distribution de la variable cible :


Fig. 1 Distribution de la variable cible

Fig. 1 Distribution de la variable cible


Distribution des co-variables en fonction de PM10 :


par(mfrow = c(3, 3), mar = c(2, 2, 2, 0) + .1)
for (i in 2:18) {
    plot(gui[,i] ~ gui$PM10, main = names(gui)[i])
}



Division des données en données d’apprentissage et de test

Description du package Caret: Le package caret (abréviation de Classification And REgression Training) est un ensemble de fonctions qui tentent de rationaliser le processus de création de modèles prédictifs. Cette librairie sera utilisé dans le cadre de l’échantillonage des données en échantillons train et test avec la fonction createDataPartition().

suppressMessages(library(caret))


Les données d’entraînement permettent d’entraîner le modèle d’apprentissage, cela signifie que que le modèle va apprendre sur ces données. Les données test permettent de tester la performance du modèle avant de le déployer. On utilise souvent 70 à 80% pour entrainer le modèle et 20 à 30% pour le tester. Si on utilise toutes les données disponibles pour l’apprentissage du modèle, on ne pourra pas tester sa performance. Ici, le ratio sera de 30/70.

set.seed(123)
partition <- createDataPartition(y=gui_without_na$PM10, p=0.7, list=FALSE)
training <-gui_without_na[partition,]
testing <- gui_without_na[-partition,]


A. Construction d’une forêt aléatoire

RandomForest est un algorithme de classification composé de nombreux arbres de décision. Elle utilise l’échantillonnage et le caractère aléatoire des caractéristiques lors de la construction de chaque arbre individuel pour créer une forêt d’arbres non corrélés dont la prédiction par comité est plus précise que celle d’un arbre individuel.


Description du package randomForest : hehe

suppressMessages(library(randomForest))


Forêt aléatoire simple

Nous allons commencer par construire une forêt aléatoire simple, sans toucher aux hyperparamètres et établir le RMSE (Root Mean Squared Error) de référence.

La formule pour calculer le RMSE est :

\[ \text{RMSE} = \sqrt{\frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{n}} \]


set.seed(123)
default_RF <- randomForest(PM10 ~ NO+NO2+SO2+T.min+T.max+T.moy+DV.maxvv+DV.dom+VV.max+VV.moy+PL.som+HR.min+HR.max+HR.moy+PA.moy+GTrouen+GTlehavre, data = training)
print(paste("RMSE =",round(sqrt(mean((testing$PM10-predict(default_RF, testing))^2)),2)))
## [1] "RMSE = 5.95"

Nous obtenons un RMSE de référence de 5,95. Cela signifie qu’en moyenne la valeur prédite varie de plus ou moins 5,96 par rapport à la valeur réelle.


Nous allons aussi calculer la valeur du R² de référence.

La formule pour calculer R² est :

\[ R^2 = \left(\frac{\sum_{i=1}^{n} (y_i - \bar{y})^2 - \sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2}\right) \]

r_squared=cor(testing$PM10, predict(default_RF, testing))^2
print(paste("R^2 =", round(r_squared, 3)))
## [1] "R^2 = 0.644"

Nous obtenons un R² de 0.64. Ce qui signifie que le taux de précision du modèle est de 64%.


Forêt aléatoire avec les paramètres ajustés

Les paramètres classiques du modèle sont appris au cours de l’aprentissage du modèle, mais les hyperparamètre doivent être définis auparavant. Dans le cas de RandomForest, les hyperparamètres pouvant être définis sont le nombre d’arbres de décision dans la forêt (ntree) et le nombre de caractéristiques prises en compte par chaque arbre lors de la division d’un nœud (mtry).

Le réglage de ces paramètres est plutôt expérimental que théorique, nous allons essayer de nombreuses combinaisons différentes et évaluer, comparer les résultats pour trouver les réglages optimaux. Dans cette partie, nous allons chercher la meilleure valeur de ntree avec do.trace et mtry en entraînant plusieurs modèles.


Choix du nombre d’arbres

Nous allons entrainer un modèle de forêt aléatoire pour prédire PM10 en fonction des autres variables avec l’argument do.trace de randomforest. Le modèle construit 500 arbres de décision et affiche des informations sur l’avancement de la construction de l’arbre toutes les 25 itérations. Ensuite, le nombre d’arbres optimal sera choisi en trouvant l’erreur OOB minimale.


set.seed(123)
RF_DoTrace <- randomForest(PM10~., data=training, ntree = 500, do.trace=25)
##      |      Out-of-bag   |
## Tree |      MSE  %Var(y) |
##   25 |    47.46    51.28 |
##   50 |    44.35    47.92 |
##   75 |    42.22    45.62 |
##  100 |    41.68    45.03 |
##  125 |     41.4    44.73 |
##  150 |    41.21    44.53 |
##  175 |    41.16    44.47 |
##  200 |    41.24    44.56 |
##  225 |    41.14    44.45 |
##  250 |    41.03    44.33 |
##  275 |    40.87    44.16 |
##  300 |    40.78    44.06 |
##  325 |    40.72    44.00 |
##  350 |    40.81    44.10 |
##  375 |    41.05    44.36 |
##  400 |     41.1    44.41 |
##  425 |    41.09    44.40 |
##  450 |    41.02    44.33 |
##  475 |    41.04    44.34 |
##  500 |    41.01    44.31 |


On remarque que l’erreur diminue de 25 à 325 arbres puis augmente de 350 à 500 arbres. D’après cette sortie, le meilleur nombre d’arbre est de 325.


Choix du nombre de noeuds


Nous allons entrainter un modèle de forêt aléatoire (ntree=325) pour chaque valeur de mtry entre 1 et 17 et calculer l’erreur OOB. Puis, nous allons représenter les erreurs sous forme de nuage de points.

nbvars <- 1:(ncol(gui) - 1)
oobsMtry = c()
for (nbv in nbvars){
  set.seed(123)
  RF <- randomForest(PM10~., training, ntree = 325, mtry = nbv)
  oobsMtry = c(oobsMtry,sqrt(mean((testing$PM10-predict(RF, testing))^2)))
}
plot(nbvars, oobsMtry)


La valeur optimale de mtry est 12 car elle conduit à l’erreur OOB minimale.


On refait un modèle avec les bons hyperparamètres et on estime ses performances.

set.seed(123)
tuned_RF <- randomForest(PM10 ~ ., ntree = 325, mtry = 12, data = training, importance = TRUE)
print(paste("RMSE =",round(sqrt(mean((testing$PM10-predict(tuned_RF, testing))^2)),2)))
## [1] "RMSE = 5.92"
r_squared=cor(testing$PM10, predict(RF_DoTrace, testing))^2
print(paste("R^2 =", round(r_squared, 3)))
## [1] "R^2 = 0.644"


Conclusion :

En ajustant les hyperparamètres, nous obtenons un RMSE légèrement inférieur (5,92 < 5.96) mais un R^2 aussi légèrement inférieur (0,642 < 0.644). Les deux modèles montrent des performances très similaires, autant garder le modèle le plus simple.

Le modèle final qui sera utilisé dans la suite du document est le modèle sans fine tuning.


Importance des variables avec randomForest

varImpPlot(default_RF,  sort  =T,  main  =  'RF  - Variable importance ')


Les variables les plus importantes sont NO2, NO, GTrouen, VV.moy, PA.moy, GTlehavre et SO2.


C. Interprétabilité avec le package randomForestExplainer


Description du package randomForestExplainer :

  • fonctions qui calculent un ensemble de mesures d’importance variable, y compris celles déjà utilisées et quelques nouvelles telles que le nombre de fois qu’un prédicteur a été utilisé pour diviser le nœud racine d’un arbre dans la forêt.

  • grande variété de possibilités de visualisation de la forêt qui aident à comprendre sa structure et à évaluer l’importance des variables.

  • fonction wrapper qui utilise toutes les capacités du package pour une forêt aléatoire donnée et produit un rapport HTML résumant les résultats.

  • une fonction permettant de sortir un fichier html avec tout les graphiques primaires du package (explain_forest()).


suppressMessages(library(randomForestExplainer))


Minimal depth distribution


La fonction min_depth_distribution permet de calculer la profondeur des prédicteurs pour chaque arbre et rassembler les résultats pour n arbres dans un seul dataframe et calcule le minimum de la profonduer de chaque variable dans chaque arbre.

#min_depth_frame = min_depth_distribution(default_RF)
#save(min_depth_frame, file = "min_depth_frame.rda")
load("min_depth_frame.rda")
knitr::kable(head(min_depth_frame, n=10), caption = "Distribution de minimal_depth")
Distribution de minimal_depth
tree variable minimal_depth
1 DV.dom 2
1 DV.maxvv 6
1 GTlehavre 5
1 GTrouen 0
1 HR.max 5
1 HR.min 7
1 HR.moy 5
1 NO 1
1 NO2 3
1 PA.moy 4


La fonction plot_min_depth_distribution utilise les informations sur la profondeur minimale de chaque variable dans chaque arbre de la forêt pour représenter graphiquement la distribution discrète de la profondeur minimale.


plot_min_depth_distribution(min_depth_frame)


Les varibles incluses dans le graphiqud sont NO, NO2, GTrouen, PA.moy,VV.moy, GTlehavre, SO2, T.moy, T.max et T.min. Les variables NO et NO2 (1.95 et 1.97) ont des profondeurs mimales moyennes plus basses, ce qui signifie qu’elles sont plus influentes dans la prédiction du modèle. En revanche, es variables T.min, T.max et T.moy (3.67, 3.79 et 3.79) ont des profondeurs minimales moyennes plus élevées, ce qui signifie qu’elles sont moins influentes que les autres, dans ces 10 variables.

En somme, ce graphique permet de faire de l’analyse de données et de sélectionner des variables dans le cadre de la modélisation des données.


Variable importance


La fonction plot_multi_way_importance utilise les données générés (datafrme contenate plusieurs mesures d’importance de variables) par la fonction measure_importance et trace deux ou trois mesures d’importances de variables les unes contre les autres.


#importance_frame = measure_importance(default_RF)
#save(importance_frame, file = "importance_frame")
load("importance_frame")
plot_multi_way_importance(importance_frame, size_measure = "no_of_nodes")


Les points sont allignés le long d’une diagonale, ce qui suggère une corrélation positive entre le nombre de noeuds et la profondeur minimale moyenne. On identifie les variables les plus importantes dans la zone supérieur droite du graphique : NO et NO2. Ce qui correspond bien avec graphique “minima depth distribution”.


plot_multi_way_importance(importance_frame, size_measure = "p_value")
## Warning: Using alpha for a discrete variable is not advised.



La fonction plot_importance_ggpairs utilise les données générés (dataframe contenant plusieurs mesures d’importance de variables) par la fonction measure_importance et représente graphiquement les corrélations entre les différentes mesures d’importance des variables.


plot_importance_ggpairs(importance_frame)


Les mesures de node_putity_increase et mean_min_depth sont fortement corrélés (-0.933). La forte corrélation négative entre ces deux variables suggère que les nœuds améliorant la pureté sont situés à des profondeurs plus élevées dans les arbres de décision, alors que les variables moins profondes ont une influence moindre sur la pureté des nœuds.


Variable interactions


La fonction plot_min_depth_interactions prend en entrée le résultat de la fonction min_depth_interactions qui contient des informations sur la profondeur minimale moyenne conditionnelle pour chaque interaction entre variables, calculé de différentes manières. Elle trace ensuite la profondeur minimale moyenne conditionnelle pour un nombre donné des interactions les plus fréquentes, ainsi que la profondeur minimale moyenne de la deuxième variable dans chaque interaction afin de comparer. Cette représentation graphique permet donc d’identifier les relations les plus importantes et de comprendre comment les différentes variables interagissent pour influencer les prédictions du modèle.


#(vars <- important_variables(importance_frame, k = 5, measures = c("mean_min_depth", "no_of_trees")))
#interactions_frame <- min_depth_interactions(default_RF, vars)
#save(interactions_frame, file = "interactions_frame.rda")
load("interactions_frame.rda")
knitr::kable(head(interactions_frame[order(interactions_frame$occurrences, decreasing = TRUE), ]), caption = "Visualisation de intercations_frame généré par la fonction min_depth_interactions")
Visualisation de intercations_frame généré par la fonction min_depth_interactions
variable root_variable mean_min_depth occurrences interaction uncond_mean_min_depth
48 PA.moy NO2 1.863839 448 NO2:PA.moy 3.118
13 GTlehavre NO2 2.407741 446 NO2:GTlehavre 3.278
38 NO NO2 2.106402 446 NO2:NO 1.954
73 T.moy NO2 2.068085 445 NO2:T.moy 3.672
42 NO2 NO 1.951982 444 NO:NO2 1.974
43 NO2 NO2 2.105290 443 NO2:NO2 1.974


plot_min_depth_interactions(interactions_frame)


  • Les cinq interactions les plus fréquentes sont NO2:PA.moy, NO2:GTlehavre, NO2:NO, NO2:T.moy et NO:NO2.

  • Les quatres interactions les plus importantes sont NO2:HR.max, NO:HR.max, NO:DV.dom et NO2:DV.dom.

  • Seule l’interaction NO2:PA.moy se situe au même niveau que la ligne rouge, cette interaction a une profondeur minimale moyenne conditionnelle similaire à la moyenne des interactions. Toutes les autres interactions ont une profondeur minimale moyenne conditionnelle plus élevée.


Prédiction de la forêt sur une grille


Nous avons calculer les interactions, maintenant nous pouvons nous intéresser à la manière dont la prédiction de la forêt dépend des valeurs des deux variables qui constituent l’interaction.

La fonction plot_predict_interaction permet de tracer la prédiction de notre forêt sur une grille de valeurs pour les composants de chaque interaction. Nous allons appliquer cette fonction aux quatres interactions les plus importantes.


NO2:HR.max


plot_predict_interaction(default_RF, gui_without_na, "NO2", "HR.max")


Le graphique montre une teinte bleue prononcée à gauche, et des nuances plus claires et du rouge à droite. Des valeurs plus élevés de NO2, conduisent à des prédictions plus élevés de PM10.


NO2:HR.max


plot_predict_interaction(default_RF, gui_without_na, "NO", "HR.max")


Le graphique montre une teinte bleue prononcée à gauche, et des nuances rouges à droite. Des valeurs plus élevés de NO, conduisent à des prédictions plus élevés de PM10.



D. Interprétabilité avec le package pdp


Description du package pdp : Ce package a pour but de construire des courbes de dépendance partielle (c’est-à-dire des effets marginaux) à partir de différents types de modèles d’apprentissage automatique en utilisant R.


library(pdp)


Les courbes de dépendance partielle (PDP) ont en effet été introduits par Friedman en 2001, principalement dans le but d’interpréter les modèles d’apprentissage automatique complexes, en particulier ceux qui dépassent la simple régression linéaire. Le graphique analyse comment la variation du prédicteur i influence les prédictions du modèle en moyennant les autres prédicteurs. En examinant comment les prédictions du modèle changent à mesure que les valeurs de i varient, on peut comprendre l’impact de cette variable sur les résultats du modèle, facilitant ainsi l’interprétation de son importance dans la prédiction.


Nous allons construire des graphiques partiels pour les variables explicatives qui sont les plus importantes dans notre modèle. Chaque graphique illustre la relation entre une variable explicative et la prédiction du modèle pour PM10, en tenant compte des autres variables dans le modèle.


partial(default_RF, pred.var = "NO", plot = TRUE)

partial(default_RF, pred.var = "NO2", plot = TRUE)

partial(default_RF, pred.var = "SO2", plot = TRUE)

partial(default_RF, pred.var = "T.min", plot = TRUE)

partial(default_RF, pred.var = "T.max", plot = TRUE)

partial(default_RF, pred.var = "T.moy", plot = TRUE)

partial(default_RF, pred.var = "VV.moy", plot = TRUE)

partial(default_RF, pred.var = "PA.moy", plot = TRUE)

partial(default_RF, pred.var = "GTrouen", plot = TRUE)

partial(default_RF, pred.var = "GTlehavre", plot = TRUE)


La valeur prédite de PM10 augmente quand les NO et NO2 augmentent. Elle augmente aussi quand SO2 augmente. La valeur prédite augmente quand T.min diminue en dessous de -5 ou augmente au dessus de 10. Elle augmente quand T.max va en dessous de 2 et augmente entre 13 à 30. Elle augmente quand T.moy est en dessous de 0 et augmente entre 10 et 25. Elle décroit strictement qd vv.moy est entre 1 et 2, puis stagne entre 2 et 10 à une valeur basse.

  • Les variables NO et NO2 ont un effet positif sur la valeur prédite de PM10, ce qui signifie que des niveaux plus élevés de ces polluants sont associés à des valeurs prédites plus élevées de PM10. De même, l’augmentation de SO2 est également associée à une augmentation de la valeur prédite de PM10.

  • Pour les variables liées à la température (T.min, T.max, T.moy), des tendances complexes sont observées : la valeur prédite de PM10 augmente lorsque T.min est inférieur à -5 ou supérieur à 10, diminue lorsque T.max est inférieur à 2, puis augmente entre 13 et 30. De plus, la valeur prédite augmente lorsque T.moy est inférieur à 0 et augmente entre 10 et 25. Quand il fait très froid ou très chaud, le taux de PM10 prédit augmente.

  • En ce qui concerne la variable vv.moy, la valeur prédite de PM10 décroît strictement lorsque vv.moy est entre 1 et 2, puis stagne entre 2 et 10 à une valeur basse.


Il est aussi possible de faire une représentation tridimensionnelle qui montre la relation conjointe de deux variables explicatives avec la prédiction du modèle, grâce à plotPartial().


pd <- partial(default_RF, pred.var = c("NO", "NO2"), chull = TRUE)
plotPartial(pd, main = "représentation tridimensionnelle entre NO et NO2 avec la prédiction")


Conclusion


En utilisant le package randomForestExplainer, nous avons exploré la distribution des profondeurs minimales des variables et évalué leur importance à l’aide de divers graphiques. Les variables les plus influentes identifiées incluent la Concentration moyenne journalière de PM10, NO, NO2, SO2, la Température minimale, maximale et moyenne journalières, la Vitesse maximale et moyenne du vent, les Pluies journalières, l’Humidité relative minimale, maximale et moyenne, la Pression atmosphérique journalière, et le Gradient de température journalier pour les stations GTrouen et GTlehavre. En examinant les interactions entre les variables avec plot_min_depth_interactions, nous avons visualisé comment elles affectent les prédictions du modèle. Enfin, en utilisant plot_predict_interaction, nous avons exploré l’impact des niveaux de polluants sur les prédictions de qualité de l’air.

En utilisant le package pdp, nous avons utilisé des courbes de dépendance partielle (PDP) pour interpréter un modèle de forêt aléatoire appliqué à des données environnementales. Les PDP permettent de visualiser comment la variation d’une variable explicative influence les prédictions du modèle tout en tenant compte des autres variables. Nous avons construit des graphiques partiels pour les variables les plus importantes telles que la Concentration moyenne journalière de NO, NO2, SO2, les Températures minimale, maximale et moyenne journalières, la Vitesse maximale et moyenne du vent, la Pression atmosphérique journalière, et le Gradient de température journalier pour les stations GTrouen et GTlehavre.